Pair-distribution functions of the two-dimensional electron gas 



o 
o 



Paola Gori-Giorgi,* Saverio Moroni, and Giovanni B. Bachelet 
INFM Center for Statistical Mechanics and Complexity and Dipartimento di Fisica, 
Universita di Roma "La Sapienza" , Piazzale A. Moro 2, 00185 Rome, Italy 
(Dated: February 2, 2008) 

Based on its known exact properties and a new set of extensive fixed-node reptation quantum 
Monte Carlo simulations (both with and without backflow correlations, which in this case turn out 
to yield negligible improvements), we propose a new analytical representation of (i) the spin-summed 
pair-distribution function and (ii) the spin-resolved potential energy of the ideal two-dimensional 
interacting electron gas for a wide range of electron densities and spin polarization, plus (iii) the 
spin-resolved pair-distribution function of the unpolarized gas. These formulae provide an accurate 
reference for quantities previously not available in analytic form, and may be relevant to semicon- 
ductor heterostructures, metal-insulator transitions and quantum dots both directly, in terms of 
phase diagram and spin susceptibility, and indirectly, as key ingredients for the construction of new 
two-dimensional spin density functionals, beyond the local approximation. 
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I. INTRODUCTION AND MAIN RESULTS 

The two-dimensional electron gas (2DEG), realized in 
semiconductor heterostructures, has been a source of 
lasting inspiration for at least two generations of funda- 
mental and applied researchers i In recent years, for ex- 
ample, a new gold rush has been triggered by the experi- 
mental discovery of a metallic phase at low temperature, 2 
in contrast with the scaling theory of localization in 
two dimensions (2D)j£ and, independently, by the scien- 
tific and technological progress on quantum dots, which, 
at semiconductor interfaces, become nothing but tiny, 
quasi-two-dimensional quantum disks, 4 

In this context, accurate predictions obtained from a 
simplified model, such as the ideal 2DEG (strictly 2D 
electrons interacting via a l/r potential within a uni- 
form, rigid neutralizing background), represent a valu- 
able reference. For example, a recent analytic represen- 
tation of quantum Monte Carlo correlation energies^ as 
a function of spin polarization £ and coupling parameter 
r s = l/^/Tmas (where n is the density and as is the Bohr 
radius) has been immediately picked by several authors, 
either because of its relevance to the phase diagram of 
the 2DEGf£ or because of the corresponding prediction 
for the spin suceptibilityf££ or, last but not least, because 
the analytic representation of the correlation energy ver- 
sus n and £ is a key ingredient for the density functional 
theory of quantum dots4i2*ifl 

Such an interest encouraged us to extend our previous 
work on energies to the spin-resolved pair-distribution 
functions g aa i (r, r') of the 2DEG, whose accuracy and 
availability in analytic form may serve a variety of pur- 
poses: the exchange-correlation hole and its dependence 
on the electron density and spin polarization may be rel- 
evant to the physics of the metal-insulator bifurcation in 
2Dii and to self-energy theories of the 2DEGji£ but is 



also needed for the estimate of the effects of the finite 
thickness on the spin susceptibility 1 ^ and for the con- 
struction of generalized-gradient approximations (GGA) 
or weighted-density approximations (WDA) of density 
functionals, in analogy to the 3D casei^i^ The availabil- 
ity of density functionals better than local-spin-density 
(LSD) approximations for the 2DEG would, in turn, al- 
low an almost exact description of quantum dots, since 
the spatial variation of their carrier density is rather 
weakAiS 

In this paper, we exploit the known exact properties 
of the pair-distribution functions (recalled in Sec. ITTft , 
and, based on a new set of extensive fixed-node quan- 
tum Monte Carlo simulations (described in Sec. IIII|) . we 
propose, in Sec. IIVI our analytic representation of (A) 
the spin-summed pair-distribution function of the ideal 
two-dimensional interacting electron gas for a wide range 
of electron densities and spin polarization, and (B) the 
spin-resolved pair-distribution function of the unpolar- 
ized gas. In Sec. [V] we discuss the quality of such an 
interpolation, and finally, in Sec. IVII , we evaluate the 
spin-resolved potential energy, of interest in the construc- 
tion of dynamical exchange-correlation potentials in the 
spin channel^ii and propose the corresponding analytic 
representation. 

As a result, quantities which are relevant to the physics 
of semiconductor heterostructures and quantum dots, 
and/or represent a key ingredient for the construction 
of two-dimensional spin density functionals beyond the 
local approximation, are now available, for the first time, 
in analytic form. 

Fortran subroutines for the evaluation of the 
parametrized quantities are available upon request to 
gp.giorgi@caspur.it or Giovanni.Bachelet@romal.infn.it. 



II. DEFINITIONS AND EXACT PROPERTIES 
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For an electronic system, the pair-distribution func- 
tions g aa i (r, r'), if n CT (r) is the density of electrons with 
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spin a 



or I, are defined as 

$|^Kr)^(rOVv(r')^(r)|$) 



'(r,r') 



7v(r)n CT '(r') 



(1) 



where ip^ and ip a are the creation and annihilation field 
operators, respectively, and is the ground-state wave- 
function. The functions g aa -i are thus related to the prob- 
ability of finding two electrons of prescribed spin orien- 
tations at positions r and r'. The normalization is such 
that the case of completely independent particles (with- 
out exchange and correlation) corresponds to the condi- 
tion g aa i = 1 . Hartree atomic units are used throughout 
this work. 

For a two-dimensional uniform electron gas, the func- 
tions g aa i only depend on r = |r — r'|, and parametri- 
cally on the density parameter r s = l/^nn and on the 
spin-polarization parameter £ = (rtf — n^jn. The total 
(spin-summed) pair-distribution function is defined as 
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For small r, when two electrons get closer and closer, 
the behavior of g aa i is governed by the cusp conditions, 18 
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r=0 



= 2g n (r = 0,r s ,0 (3) 

= g aa {r = 0,r„C) =0 (4) 
d 2 

= 2 Q^9aa{r,r S7 _ ■ (5) 



S acr ,{q,r s ,() = <W + 



dx [g aa > - 1] x J (qx), 
(6) 

where q = k/k? is a scaled variable in reciprocal space, 
and Jo is the Bessel function of order 0. The total (spin- 
summed) static structure factor is 



S 
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(7) 



its long- wavelength (i.e., small-q) behavior is determined 
by the plasma collective mode^ 



S(q^0,r s ,C) 
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0(q 2 ), 



(8) 



and thus does not depend on C- 



Usually g aa i (and consequently S aa i ) is conventionally 
divided into the (known) exchange and the (unknown) 
correlation terms, 
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(9) 
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where J\ is the first-order Bessel function, 9 is the Heav- 
iside step function, and k\ — fep-y/l + (, k F = kp\/Y~~(. 
The functions g x and S x correspond to a uniform 2- 
dimensional system of noninteracting fermions; once the 
scaled variables x and q are used, they do not depend 
explicitly on r s : g x — g x (x,<^), S x = S x (q,()- In what 
follows, we use the name pair- distribution function for 
the whole thing (g = g x + g c , exchange plus correlation), 
and pair- correlation function for its correlation-only con- 
tribution g c . 

Combining Eqs. JSJ) and 1141) . we find the small-g be- 
havior of the spin-summed correlation static structure 
factor, 



Eqs. (J3J) and (f5j> are due to the dominance of the potential 
term l/|r — r'| in the many-body hamiltonian as r — > r'; where 
Eq. J3J comes from the Pauli principle. 

At this point, it is convenient to introduce the scaled 
variable x = kpr, where kp = \2/r s is the Fermi 
wavevector of the unpolarized gas. 

The Fourier transforms of g aa i — 1 are the spin-resolved 
static structure factors*^ which, for a 2D uniform gas, 
are 



S c (g^0,r fl ,C) = — 0(C) 9 + 

7T 



3/2 

+ 0(q 2 ), (15) 
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plays the same role of the three-dimensional function (f> 
of Refs .12(1 and l2ll As well known from the properties of 
Fourier transforms, the small-g behavior of S determines 
the oscillation-averaged long-range part of g. We thus see 
that, individually taken, g c and g x — 1 have long-range 
tails oc r -3 ; but these tails exactly cancel in the pair- 
distribution function (exchange plus correlation), so that 
g — 1 = g x + g c — 1 approches zero as r~ 7 / 2 . 

While the long- wavelength limit of the total S, Eq. |(HJl, 
is well known, little is known about the small-g behavior 
of the spin-resolved S aiJ i (and hence about the long-range 
part of gcrcr' )■ The conservation of the number of particles 
implies 



Saa'(q = 0,r s ,C) = 0. 



(17) 



In Section flV Bl we discuss an approximate expression for 
S aa '(q — > 0, r s , C = 0) consistent with our QMC results. 

Finally, the spin-summed g c yields the correlation part 
of the expectation value of the Coulomb potential energy, 
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v c (r s , £), which can be obtained from the correlation en- 
ergy e c (r s , C) via the virial theorem^ 

y / dxg c (x,r s ,C) = v c {r s ,() = — 7^7 [rf e c (r s , 0] • 

(18) 

III. QUANTUM MONTE CARLO 
CALCULATION 

The ground-state expectation value (5 of a local opera- 
tor O, such as the pair-distribution function or the static 
structure factor, is estimated as 

= (#(j3)\6\*tf))/(*(P)\*(J3)) (19) 

using a reptation quantum Monte Carlo (RQMC) 
algorithmic Here \1/ is a trial function, and = 
c -PH/2^ can mac j e sufficiently close to the exact 
ground state $ by choosing the "imaginary time" /3 large 
enough. 

The estimate of Eq. ljT§|) is called "pure" , as opposed 
to the "mixed" estimate O mix = usually 
adopted in connection with the diffusion Monte Carlo 
(DMC) method, 24 More precisely, previous DMC re- 
sults for the pair-distribution function of the 2D electron 
gaa25i2Si2I have been based on extrapolated estimates, 
6 ext = 26 mlx - (tf 1 6 1 *}/(tf|tf>. The bias in 6 ext is 
quadratic in the error of the trial function. Such an es- 
timate is often very accurate, but a well converged pure 
estimate, as obtained in the present work, has the ad- 
vantage of being independent of the quality of the trial 
function \? (except for its nodal structure, see below). 

The RQMC method features a discretized path inte- 
gral representation of the importance-sampled imaginary 
time propagator, 

G(R -» R P ;P) = y(Rp)(R P \c-P H \R )/y{R Q ) (20) 

— J dRi ■ ■ ■ dRp^i!lf =a 1 G(Ri — > e), 

where e = (3/P is the time step and Ri is the set of the 
2N coordinates of the N electrons at the i-th step. We 
use the standard short-time approximation 2 

G(R^R';e) ~ Ae -[R'-R-evmnm 2 /2e (21) 

x e -e[E L (R')+E L (R)]/2^ 

where E L (R) = H^(R)/^(R) is the "local energy", and 
A = (2Tte)~ N is a normalization constant. Replacement 
of Eqs. (|21|l and 11' 21 into (|19|l yields an integral amenable 
to Monte Carlo evaluation, using a generalized Metropo- 
lis algorithm to sample paths in an enlarged configuration 
space, X — {i? , ■ • • , Rp}- 

In our simulations wc consider iVf spin-up and N± 
spin-down particles in a square box with periodic bound- 
ary conditions. The spin-resolved pair-distribution func- 
tions are obtained 2 ^ averaging VdN aa >(r')/[N a (N a i — 



d aa ')2iTrA] in the middle slice of the path during the 
simulation, where V is the volume of the simulation 
cell, dN aa i{r') is the number of electron pairs with dis- 
tance r' between r — A/2 and r + A/2. The struc- 
ture factors are computed analogously, for vectors k in 
the reciprocal lattice of the simulation cell, by averaging 
p CT (k) jOCT ,(-k)/(7V CT 7V CT 1/2 , where p a (k) = £\ exp(-*k • 
Yj) is the density fluctuation of electrons with spin a. 
The total number of particles is 42, 50, 50 and 45 for po- 
larization 0, 0.48, 0.80 and 1, respectively. By repeating 
simulations for different system sizes in the unpolarizcd 
case, finite size effects on the pair-distribution function 
have been estimated to be of order 0.01. The systematic 
bias due to finite projection time and finite time step can 
be kept within this level by suitable choices of the param- 
eters (5 and e. In our simulations, this results in paths of 
501 time slices. 

We avoid the fermion sign problem using the fixed node 
approximation (FNA)^ whereby the paths are not al- 
lowed to cross the nodes of the trial function. The FNA, 
which gives the lowest energy upper bound consistent 
with the nodal structure of the trial function, is the only 
source of uncontrolled approximation in the present cal- 
culation. In order to gauge the sensitivity of the com- 
puted pair-distribution function on the nodal structure 
of 5', we have performed our simulations using two trial 
functions with different nodes. 

Our first trial function is of the simplest Jastrow-Slater 
form ^{R) = J(R)S(R). Here J(R) =U l<3 : exp(-u(r tf )), 
Tij being the distance between the electrons i and j, is a 
symmetric Jastrow factor; it describes pair correlations 
through the function u(r), which is optimized (by mini- 
mizing the variational energy) for each density and polar- 
ization; it's always positive, so it does not alter the nodal 
structure, which is entirely determined by the other fac- 
tor S(R), a product of two Slater determinants (one for 
each spin component) of plane-wave one-particle orbitals 
exp(ki -Tj). 

Our second trial function has the same Jastrow factor, 
but its nodal structure is more accurate, since it includes 
"backflow" correlations 2 ^* 2 ^ by replacing the electron co- 
ordinates Tj in the Slater determinants with "quasi- 
coordinates" 

X J = r 3+J2 Vi^i ) ( r * - r i)> ( 22 ) 

where rj(r) is another function to be optimized for each 
density and polarization. 

In a previous variational calculation 2 ^ the difference 
between the pair-distribution function calculated with 
the simple Slater- Jastrow and the backflow trial func- 
tion was found to be of order 0.03. Here we find that, 
in a fixed-node calculation, such effect is even smaller: 
Figure n shows the difference in and g^, computed 
with cither plane- wave or backflow nodes, for r s = 2 and 
20 at zero polarization. In the worst case (large r s , lower 
panel) these differences are half as large as found in the 
variational casej2I while for small r s (upper panel) they 
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FIG. 1: Numerical difference between the spin-resolved pair- 



distribution functions g 



PW 

9 (7<t' 



at r s 



2 (upper panel) 



and r a — 20 (lower panel), as obtained from two fixed- node 
simulations with different nodal structures. The superscript 
indicates backffow (BF) or plane-wave (PW) nodes. 



are much smaller than that. These differences are essen- 
tially invisible on the scale of our g aa i (r, r s , Q calculated 
with plane- wave nodes, some samples of which are shown 
Figure|21 As a consequence, an analytic representation of 
the spin-summed g(r, r s , £) and of g a cr'( r , r s,C — 0) (see 
next Sec. HVf) based on the plane-wave results, as the one 
presented here, happens to give an equally good repre- 
sentation of the backflow results, because the difference 
due to the improved nodal structure is either comparable 
or smaller than the fitting error. 



IV. ANALYTIC REPRESENTATION 

In this section we describe our analytic representations 
of the spin-summed pair-correlation function g c (x,r s ,() 
valid for 1 < r s < 40 and < C < 1> an d of the spin- 
resolved g c aa , (x, r s , () for £ — and 1 < r s < 10. These 
functions arc built along the lines of R.efs. 20.21 and of 
Ref. for the 3D case. 

The strategy is the following. We build the spin- 
summed g c as a sum of three terms: long-range, short- 
range, and oscillatory. The long-range term is taken from 
the random-phase approximation (RPA) and multiplied 
by a cutoff function which quenches its short-range con- 
tribution. The short-range part is built according to the 
cusp conditions of Eqs. ©-©, as a weighted sum of ff, 
11 and || terms which, in turn, have been determined 
for £ = by a fitting procedure to the QMC results. 
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x = k c r 



x = k c r 



FIG. 2: Sample of spin-resolved pair-distribution functions as 
directly obtained from our QMC simulations (no fitting here). 



For C 7^ 0, an exchange- like ^-dependence of these a a' 
short-range coefficients has been assumed. The oscilla- 
tory part is empirical, being entirely determined by a 
fit to the QMC data. The analytic function g c is also 
constrained, via Eq. \V&\. to reproduce our parametrized 
correlation energy of Ref. 0. 

The analytic parametrization of the spin-resolved g c aa , 
is more difficult, because less is known about its exact 
properties. We had to rely more heavily on our QMC 
data, and, for the time being, we successfully interpo- 
lated <7|| only in the unpolarized case = 0) and for 
r s G [1,10]. This parametrization, combined with the 
one for the total g c also yields gS^ — gf* — 2g c — g^, . We 
build <?Si using a functional form similar to the one just 
described for the total g c : a sum of a long-range term, 
a short-range term, and an oscillatory term. The long- 
range term is obtained by a modification, consistent with 
our QMC data, of the long-range analytic form appropri- 
ate for the total (spin-summed) g c . The short-range term 
is simply the || part of the total g c . The oscillatory part 
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FIG. 3: Spin-summed pair-distribution function (exchange 
plus correlation: g = g x + g c , see text) for four different 
values of the spin-polarization parameter £, and for r s = 1, 2, 
5, 10, 20, and 40 (larger r a values have stronger oscillations). 
The dots correspond to our QMC data, the solid lines to our 
analytic representation. Error bars are comparable with the 
dot size. 



is, again, empirical. Because the short-range parts of the 
total g c and of share some parameters, we performed 
a simultaneous, global, three-dimensional (x,r s ,() fit of 
g c (x,r s ,() and g^(x,r s ,Q — 0). This procedure and all 
the relevant equations are detailed in the next subsec- 
tions. 



A. Spin-summed pair-correlation function 

We parametrize the spin-summed g c as 

6 

9 C = [9lr{x) + g OS cm(x)} F cut (x) + ^ c n x n , (23) 



where <7lr is a long-ranged function whose Fourier trans- 
form exactly recovers Eq. I|15l) . .g sciii is an oscillating 
function to be fitted to the QMC data, and the last term 
in the r.h.s takes care of the short-range properties. The 
function _F cut (x) quencheaSi the short-range contribution 

Of (5LR + ffoscill), 



F cut = 1 - e- dx ~ (1 + dx A + t^x 4 + U'x' 



2 i 1j2 4 , 1j3„,6^ 



(24) 



The parameter d(r s ) determines the mixing of long-range 
and short-range terms in Eq. (|23[) . 



1. Long-Range part 

The long-range part is built with the same procedure 
used for the 3D case in Refs. I20II21I and l38l and detailed 
in Appendix lAl 



(25) 



where v = ^/2r s (f> 2 x is another scaled variable, and 4>(Q 
is given by Eq. (|16|) . The function fi(v) is reported in 
Appendix ^ 



2. Short-Range part 

The short-range part of our g c is the last term in the 
r.h.s of Eq. $2$. We have 
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(26) 
(27) 
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(28) 



c 3 = da 
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where a™ are the short-range coefficients of the spin- 
resolved pair-distribution functions, 



W (z^0,r s ,C)=5>r>, 



and 



3&F 



(30) 



(31) 
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The pair-correlation function at zero electron-electron 



distance, or "on-top" value 5^(0) 
parametrized as 



,T1 



5^(0) = [1 + (a- 1.372) r, + b r 2 + cr 3 s ] e 



1, has been 



1. (32) 



The parameters a = 1.46, b = 0.258, c = 0.00037 are 
fitted to the QMC results; the exact high-density slope, 
1.372, is taken from Ref. I39L 

As said at the beginning of this section, we determine 
the spin-resolved short-range coefficients for the ( = 
case, and then we assume an exchange- like ^-dependence. 
This means that in Eqs. l(2T)|) - (|2l?|> the values of g^(0), 



,11 



and only depend on r s (not on £), and that the 



coefficients aJ^ and have the simple £ dependence 



,11 



a| T (r s ,C) = i( 1 + CK(^) 5 



(33) 



,11 



(r s ,C)=4 T (r s ,-C) 



with 

The linear parameters C4 and C5 will be used to 
constrain g c to yield the correlation energy of Ref. 
and to fulfill the particle-conservation sum rule [S(q = 
Q, r sX) = OL as m Ref- l2ll The parameter c 6 (r s ,Q is 
used to give more variational freedom to our g c for an 
accurate fit of the QMC data at higher r s . 



3. Oscillatory part 

The oscillatory part of our g c is similar to the form 
used by Tanatar and Ceperley^ 



ffosciii = — — r e _r ™ 2 x cos(m 3 x + m 4 ) (34) 
x + 1 



which is able to accurately fit the QMC data at low densi- 
ties. The exponential cutoff ensures that (fosciii does not 
alter the long-range properties embedded in (?lr- The 
parameters rrii depend on both r s and £. 



4- Sum rules 

As said, the role of the parameters C4 and C5 which 
appear in the short-range part of our g c is to fulfill the 
normalization sum rule [S c (q — 0) = 0] and to recover 
the correlation energy e c (r s , £) of Ref. yj We obtain 
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(35) 
(36) 




FIG. 4: Spin-summed pair-correlation functions g c from our 
analytic representation. Upper panel: for £ = 0, we show 
g c for r s = 2, 3, 4, 5, 7, 10, 15, 20, 30, 40; stronger oscillations 
correspond to higher r s values. The solid lines correspond 
to r s — 2, 5, 10, 20, 40, for which our g c accurately fits the 
QMC data; the dashed lines are the results for intermediate 
values of r s . Lower panel: for r s = 2, we show g c for different 
values of the spin-polarization £ = 0,0.3,0.48,0.7,0.8,0.9, 1; 
more negative "on-top" values g c (x = 0) correspond to lower 
values of £. The solid lines correspond to ( — 0,0.48,0.8, 1, 
for which our g c accurately fits the QMC data. The dashed 
lines correspond to intermediate values of £. 



with 



SLR 

Soscill 
-^oscill 



£2 30F 3 c 6 

2d 4tfV 2 2d 2 ~ ° 3 
+2<j) 5 r 2 s LR - s oscin 



C3 



5/2 d i 

(37) 
15 v 7 ^ 



2y/d 2d 4tfV 2 2d 2 
-2cj) 5 r 2 E LR - E oscm + V2r s v c 

fi(v) [1- F cut (x)]dx 



C6 



fi(v) 



F cut (x) dx 



g osc m(x) x F cut (x) dx 



g OS cm{x)F cut (x)dx, 



16 d 7 / 2 
(38) 

(39) 
(40) 
(41) 
(42) 



and v c {r s ,Q given in Eq. JTBJ. Equations l|3Ti )l -l|32 )) are 
evaluated numerically for given r s and 
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= 


-0.0151 


(3) 
Pi 


= 2.14 


(3) 
?i 


= 0.394 




= 0.045 


(3) 
?2 = 


-0.0299 


(4) 
Pi 


= 6.39 


(4) 

?i 


= -0.592 




= 2.7-10" 


-4 (4) _ 
?2 — 


-1.8 • 10" 4 



9m- 



J? 

(2) 
(3) 

„( 4 ) 



: 1.1 

0.479 
0.6 
1.99 
: 1.437 



7^ =29 

,(!) 



4 3) 

,( 4 ) 



0.1 



TABLE I: Optimal parameters for the analytic representation 
of g c (x , r s , Q and g? , (x, r s , £ = 0) as described in Sec. II VI 



5. Fitting parameters 



The parameters d(r s ) , 



(r s ), al l (r s ), a p (r s ), 



C6(r s ,(), mi(r s ,C) are used to fit the QMC data. Their 
r s and £ dependence is smooth and allows for an analytic 
representation of g c (x,r s ,() valid at all r s £ [1,40] and 

Ce[0,l]: 



d(r s ) 

4 l (rs) 
alHr s ) 
a p (r s ) 

ce(r s ,C) 
m 1 (r s ,C) 

m 2 {r s ,0 



m 3 (r s X) 
m 4 (r s ,C) 



1 + 6 2 r% 

(1-X irs +X 2 r 2 s )e- X ^ 

7 (6) (C)e -7< 6) (0/^ 

/4 2) (C) 

l + /4 2) (C)r, 

M i 3) (C) + 2.7 M f (CK 



l + ^ 3) (C)r 



^ 4) (C)+5.36^(Qr 



(4), 



i + /4 4) (CK 2 



(43) 

(44) 
(45) 
(46) 

(47) 
(48) 

(49) 
(50) 
(51) 



The functional form of the short-range coefficients aZ" 
is very similar to the one used for the 3D case in Ref. l40t 
the corresponding parameters are determined by simul- 
taneously fitting the data for g^ (see next Section) and 
those for the total g c . The parameter cq only comes into 
play at high r s : its functional form, Eq. I|47|l. makes it 
vanish very rapidly as r s decreases. The same argument 
applies to the oscillatory part, whose magnitude is de- 
termined by the parameter mi of Eq. H48|). The low- 




FIG. 5: IJ. pair-distribution function (exchange plus correla- 
tion, see text) for £ = and r s = 1, 1.5, 2, 3, 5, 7, 10 (the 
larger r s values have stronger oscillations). The dots corre- 
spond to our QMC data for r 3 = I, 2, 5, 10; the solid lines 
is our analytic representation at the same r s -values. Dashed 
lines correspond to our analytic representation for the other 
values of r s . 



density limit of the parameters 7713 and r/14, 2.7 and 5.36 
in Eqs. (|50|l and (|51|l . are taken from an oversimplified 
model of localization on the sites of a triangular latticed 



The C dependence of the parameters 7, 
represented by a quadratic form: 



(6) 



7, (6) (C) = Pi + mC 



In) . (n) yx 

Pi' + <U ' C ■ 



and /if 1 ^ is well 

(52) 
(53) 



The final 32 free parameters (plus 9 parameters for gS,, 
detailed in the next section) are fitted to our data set 
(100 values of x for each r s = 1,2,5,10,20,40 and 
C = 0,0.48,0.8,1 plus those for g^ at C = and 
r s = 1,2,5,10 - a total of 2800 data), and are reported 
in Table [Q 



B. Spin-resolved pair-correlation functions = 0) 

We parametrize the fj, correlation function with a 
functional form similar to the one used for the spin- 
summed g c , 

5 

9 C U = lAx)+9lL m (x)} F cut (x)+e- d * 2 ]T <#x n , (54) 



where the function F cxli (x) and the parameter d(r s ) are 
given in Eqs. and |@3J|, respectively. 



1. Long-Range part 

While the long-range part of the spin-summed g c , 
Eq. (|25p. can be obtained from RPA, the spin- resolution 
is more problematic. Nonetheless, RPA can give some 
hints j2i especially in the r s — > limit. From RPA we 
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obtain, up to 0(q 2 ), 



3. Oscillatory part 



(55) 



with Cn(0 = l, Ctt(C) - 2/VTTC - VW/VT+C, 

and Cu(0=^n(-0- 

Here, we only treat the ( — case, for which we 
also produced spin-resolved static structure factors with 
QMC. We write the small-g part of S° a , as the RPA re- 
sult plus an r s -dependent correction, similar to the 3D 
case^ i.e., up to 0(q 2 ), 



For <7| , we use the same form [Eq. |2] of the total g c , 

: cos(ml l x + m\ l ). (63) 



,T1 



X+l 



The parameters ml depend on r s and are fitted to the 
QMC data. 



4- Sum rule 



S° tr ,(q^0,r s ,t = 0) = -q 



1 / \ 

- + a acr ,{r s ) 

TT 



7 3/2 



(56) 

with a-\--\(r s ) = — a||(r s ). This small-g behavior embod- 
ies the following properties: (i) the corresponding spin- 
resolved pair-distribution function g aa i (r) are more long- 
rangeo™ than the spin-summed g(r), (ii) parallel- and 
antiparallel-spin correlations give identical contributions 
to the plasma collective mode. The correction a^i(r s ) 
has been determined from the QMC results in reciprocal 
space for 1 < r s < 10, and is well represented by 



a n (r s ) = 0.00914 r s . 



(57) 



Thus, for the spin-resolved long-range part we use a scal- 
ing law similar to the one of Eq. 125f) . 



the function fi(v,a) is described in Appendix IB1 



(58) 



2. Short-Range part 

The short-range part of gi, is the | J. part of the total 
g c [see Eqs. We thus have 



Jl 



J ! 



'2 

J ! 



5n(0) 

4 [5u(0) 



del 1 



,11 



dc[ 



z 3 ' 



(59) 

(60) 

(61) 
(62) 



The sum rule Eq. I|17|> determines the linear parameter 



J I 

c 4 J 



with 

s U - 

LR _ 

s U - 

oscill 



c 4 



ji Jl 



"0 



c Vt c 



U 

2 c n 



3^ 



2d 4d 3 / 2 2d 2 ~ 3 8d 5 / 2 
+2^r s 2 ^ R -^ cm , 

OO 



fi(v,a n ) [1 -F cut (»)]dx 



n 



(64) 

■n 15 \^ 

^ 5 16 d 7 / 2 
(65) 

(66) 

(67) 



5. Fitting parameters 

From the global fit described for the total g c , we also 
find the r s dependence of the coefficients Cg and m\ ■ 



m\ l (r s 
m\ l {r s 



(i) 

i + 4 X V S 

„(2) 



.(3) 



^r 2 
1 + zfV 2 



,(4) 



m^(r £ 

The values of 7^ and i/j^ are reported in Tabled 



(68) 
(69) 
(70) 
(71) 

(72) 



where fff^(O), a 2 , and 03 are given in Eqs. lO, (EH, 
and (1431) . respectively. 



1 1 

The linear parameter C4 is used to fulfill the normal- 
ization sum rule of Eq. (|17fl : the parameter cl^(r s ), in- 
stead, increases the variational flexibility of , and is 
fitted to the QMC data. 



V. RESULTS 



A pictorial evidence of the quality of our analytic rep- 
resentation clearly emerges from Fig. |31 where we show 
our analytic representation for the spin-summed g(r), to- 
gether with the corresponding QMC data, for r s = 1, 2, 
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-0.1 
-0.2 
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C = o ■ 
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q = k/k F 



FIG. 6: Upper panel: total (spin-summed) structure factor as 
directly obtained from our QMC simulations (data with error 
bars) and as a Fourier transform of our analytic representation 
of g c (solid lines), for r s = 1,2,5,10,20,40. Higher peaks 
correspond to larger r s . Lower panel, same comparison for 
the IJ. static structure factor: the data with error bars are 
QMC simulations and the solid lines are Fourier transforms 
of our analytic g?\. Here the r s values are 1, 2, 5, and 10 and 



the larger deviations from the noninteracting value Sfi 
correspond to larger r s values. 







CO 



CO 

I 



CO 




5, 10, 20 and 40 and four different values of the spin- 
polarization Figure^] instead, shows that our analytic 
g c (r, r s , C) smoothly interpolates the QMC data not only 
as a function of x — kpr, as e.g. shown in the previous 
Figure |3J but also as a function of r s (upper panel) and 
of C (lower panel) . Figure El summarizes similar results 
for <7|| at C = 0. The static structure factors for ( = 
are reported in Fig. |SJ In the upper panel, we compare 
the total S(q) corresponding to our analytic g c with our 
QMC calculation (see Sec. IHIjl : the agreement indicates 
that the long-range part (q — > limit of S) of the analytic 
g c has been accurately described. In the lower panel, we 
show similar results for S-^i(q). We see that the long- 
range (q — ► 0) spin-resolution of Eq. I|56|l is consistent 
with the QMC results. 

Recently, Atwal, Khalil and Ashcroft^i (AKA) have 
presented a parametrization of the dynamical local- 
field factors (spin-symmetric, |T + IJ-i an d spin- 
antisymmetric, IT — 41) for the C = 2D electron gas, as 
a function of the wavevector q and of the imaginary fre- 
quency iuj. Following the analysis carried on for the 3D 
electron gas by Lein, Gross and Perdew^ Asgari et al£& 
have compared the wavector decomposition of the cor- 
relation energy resulting from the AKA spin-symmetric 
local field factor with the one resulting from our present 
work, based on QMC results, for r s = 1. In the upper 
panel of our Fig. [3 we make a similar comparison (in this 
case at full coupling strength) for r s = 2 and r s = 5. 
In the lower panel of the same figure we also compare 
the results from the AKA spin-antysimmetric local-field 
factors. We see that the spin-summed AKA S c (q) is in 
fair agreement with our result for q < 1.5, where both 
curves recover the exact behavior of Eq. I|15|) . The spin- 
antisymmetric AKA curves are, instead, quite different 
from our result, even for small q. This discrepancy prob- 
ably comes from an inaccurate description of the high-cj 
behavior of the q — > limit of the AKA parametrization 
for the spin channel^ In particular, Eq. (26) of AKA 
yields a formally divergent result when combined with 
the known limiting behavior— Sjl(9 — > oo) oc q~ 3 . 



VI. SPIN-RESOLVED POTENTIAL ENERGY 

The correlation part of the potential energy, v c (r s ,£) 
of Eq. (|T%)l . can be divided into ft, II, and || contribu- 
tions, such that v c = vl* + vjr- + vl^ . These spin-resolved 
components of v c are important ingredients for the study 
and construction of dynamical exchange-correlation po- 
tentials in the spin channel^ii They can be written as 
the expectation value of the Coulomb potential 1/r on 
the spin-resolved g c a(J i , 



FIG. 7: Upper panel: spin-summed correlation static struc- 
ture factors from our analytic representation of g c and from 
the Atwal, Khalil and Ashcroft^- (AKA) dynamical local-field 
factor. Lower panel: the same comparison is done for the spin 
channel (correlation only, jj — ||)- All curves are for C, — 0. 



(r.,C) 



(2 - tW) n a n i 
V2r s n 2 



9aa' (x,r s ,()dx. 



(73) 

We have evaluated the r.h.s. of Eq. (|73|) by numeri- 
cal integration of our QMC data for g c aa , (x, r s , C), at 
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FIG. 8: Relative error Av c (r s , Q/v c (r s , <) = (w c - n™ T )/^ 
between v c calculated using the r.h.s. of Eq. I|18|l (with e c 
from Ref. 0), and i>J NT , obtained by numerical integration of 
g c (see text). 




5 10 15 20 25 30 35 40 



FIG. 9: Fraction of jj, contribution to the correlation part of 
the potential energy, F-\\ = vl 1 ' /v c . 



sented by 

F u (r s ,C) = ^ T H T D (C) + ^i(C)^+^(C)^] x 



x log 1 



MO 



where the high-density FH D is given by Seidl^ 



^T H T D (0 = 

^(0 = 



-19.54(1 + 



(75) 



(76) 



153.38^(0 - 192.46' 

(i + C)iog(i + C) + (i-C)iog(i-C) 

21og(2) 

+0.0636 C 2 - 0.1024 ( 4 + 0.0389 C 6 , (77) 

and the functions u)i(C) have been obtained by htting 
our data for w] T (r s ,C) for C = -0.8, -0.48, 0., 0.48, 0.8 
(the negative C values corresponding to the || data), 

wi(C) = (1 -C) (-0.006- 0.03 C) (78) 
MO = (1-0 (-0-01 + 0.03C) (79) 
MO = 3.6 (1 + 4 - (80) 

Equations H75|) - (|80|) completely determine the spin res- 
olution of v c (r s ,C), since F u (r s ,C) = F TT (r s ,-C) and 

F 1l = 1 - ^TT - F ll- 

In Fig. El we show our numerical results for the 
antiparallel-spin fraction F^(r s , 0, together with our fit- 
ting function; the relative errors on the fit of (not 
shown) are of the same order of magnitude of those of 
Fig. 03 We see that the correlation part of the potential 
energy is completely dominated by the fj. contribution, 
even for £ as high as 0.8. 



C = 0,0.48,0.8, and r s = 1,2,5,10,20,40. This means 
that the integration in the r.h.s. of Eq. I|73(l has been 
truncated at L/2, where L is the side of the simulation 
cell (in our case L/2 ~ 6r s ). The resulting i£ cr are thus 
affected by the finite-size error, since they correspond to 
systems with fixed number of particles (see Sec. lIIIf) and 
an infinite-size extrapolation is not available in this case. 
One can get an idea of the magnitude of such error by 
using the same numerical-integration procedure for the 
spin-summed g c , and then comparing the results with 
the corresponding thermodynamic limit, the last term of 
Eq. (|18fl . combined with our— e c (r s , Q. The relative error 
between the two evaluations of v c is reported in Fig. [S] 
it is of the order of few percents. At £ = 1, Figure 
disproofs Eq. (10) of Ref. H, which predicts a different 
virial theorem for the fully polarized system. 

We have parametrized our spin-resolved v" 17 (r s , 



t£ a '(r a ,Q=F aa ,(r a ,Qvc(r a ,Q. (74) 
The fractions F aa (r s , f° r parallel spins are well repre- 
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APPENDIX A: LONG-RANGE SCALING 

In this appendix we describe the construction of the 
long-range part of our g c . We follow the same procedure 
used for the 3D case in Ref. |38|, where the interested 
reader can find more details and comments on the rele- 
vant physics. Here we briefly recall the main equations 
and emphasize the differences between the 3D and the 2D 
case. As for the 3D case^ the results of this appendix 
can be used in the construction of a generalized-gradient 
approximation for a 2D correlation energy functional. 

Following Ref. |3^, we seek a scaling law for the 
long-range part of g c . We call "long-range" part the 
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FIG. 10: The function f(z,Q [see Eq. JA1H . evaluated within 
RPA, for three different values of the spin polarization param- 
eter f. Also shown: the exact small-z behavior of Eq. JA3II 
and the function h(z) of Eq. I|A8(I . 



oscillation-averaged asymptotic behavior of g c for large 
x = kpr, or equivalently, the behavior of S c for small 
q = k/kp. From Eq. I|15fl . we see that such a scaling law 
has the form 



S c (q^0,r s ,O^V2r s 3 /(z,C), 



where 



ZTF<, 



l V2i 



(Al) 



(A2) 



is a variable on the scale of the Thomas-Fermi wavevec- 
tor, kxF (which does not depend on r s in the 2D case), 
and the function /(z, has the small-z expansion (inde- 
pendent of () 



°>0 = - ^ + ^z 3 / 2 + 0(z 2 ). 



(A3) 



The random-phase approximation (RPA) exactly recov- 
ers Eqs. HA1 |I -HA3 [I . As in the 3D case^ 8 . we can thus 
obtain the function /(z, £) from RPA. Its (wrong) short- 
range behavior will be quenched in our parametrization 
of g c by the cutoff function F cut (a;) of Eq. 

We thus evaluated the function /(z, £) via the standard 
RPA equation 



<Srpa(<7,^C) 



du 



03 T +/3l) 2 



qk F /7r-(^+/3 l )' 



with 



(3(q,u>) = 

irq 



and 




(A4) 



(A5) 



/3 T =/3 



-c, 

(A6) 



The resulting f(z, C) has the small-z expansion of 
Eq. HA3fl . and for large z is equal to 



/(*-oo,C)=j1(C) + 



B(C) 



(A7) 



Equation (|A7|I corresponds, in real space, to a divergent 
short-range behavior. As in three dimensions, 38 with a 
suitable cutoff Eq. (|A7|) recovers the RPA high-density 
expansion of the correlation energy, e^ PA (r s — > 0,£) = 
orpa(C) + ^rpa(C)?" s log^s + 0(r s ). However, unlike the 
3D case, in two dimensions the RPA correlation energy is 
not exact even in the r s — > limit. 37 Thus, while in three 
dimensions22i2ii2 8 i the large-z behavior of /(z, £) was im- 
portant to recover the exact high-density limit of e c , there 
is no need in 2D to keep Eq. I|A7(1 in our parametriza- 
tion. Moreover, as in three dimensions^ we find that 
the C dependence of /(z,C) is very weak (see Fig. Utty . 
so that we can replace f(z,() with /(z, 0) [this is exact 
in the "important" part of /, i.e., the small-z regime of 
Eq. 1A3() ]. We thus define the function h(z) 



h(z) = f(z,0)-A(0) 



B(Q) 



(A8) 



where A(0) = -0.272076 and 5(0) = 10/tt - 3 corre- 
spond to the large-z expansion for £ = of Eq. (|A7|) . As 
shown in Fig.EH the function h(z) has the same small-z 
behavior of /(z, £), but goes to zero when z — * oo, which 
corresponds to a less diverging short-range part in real 
space. The Fourier transform of h(z) defines the function 
h{v) Of Eq. G3, 



h(v) 



h(z) z Jq(vz) dz, 



(A9) 



where v = y/2r s (j) 2 x is the appropriate real-space scaled 
variable. The function fi(v) has been evaluated numeri- 
cally, and then parametrized as 

&1 v 1/2 + b 2 v + b 3 v z l 2 + b A v 2 + be, v b l 2 



fi(v) 
with 

&6 = 
65 = 



-b 4 v 2 

{V 2 + 62)5/2 



2/tt 
9 1 



K v 
(A10) 

(All) 
(A12) 



64 = —3 60 



bi 



2b, 



-22 



1/2 



2b. 



B 



-</2 



"11 



3b 2 



2b, 



3/2 



5 5 

4' 4 



B -,- 



3 7 

4' 4 



(A13) 



61 



-64 



3.46. 



Eqs. (IA11I) - (|A13() guarantee that the Fourier transform 
of fi(v)/v satisfies Eq. (IA3|) . r(x) and B(x,y) are 
the standard Gamma and Beta functions^ 3 - T(3/4) « 
1.225416702, B(3/4,7/4) « 0.8472130848, B (5/4, 5/4) « 
0.6180248924. 
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APPENDIX B: SPIN-RESOLUTION OF THE 
LONG-RANGE PART « = 0) 

The long-range part of g° a ,(x, r s , £ = 0) has been 
simply approximated with Eq. I|58|l . where the function 



fi(v,a) is obtained from fi(v) of Eq. (|A10|) by replacing 
be with &6 = 2(l/7r + a), and consequently changing 64 
according to Eq. (|A13|I . In this way, the corresponding 
S^ a ,(q, r s , C = 0) exactly recovers Eq. 
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